{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 时间序列\n",
    "1. 平稳性Danile检验(base on Spearman秩相关系数) 418页\n",
    "2. AR自回归预测模型. 经典预测模型有: 生长曲线, 指数平滑法. 这两种对于短期波动把握不高.  \n",
    "AR自回归模型考虑了时间序列上的依存性, 又考虑了随机波动的干扰性, 对短期趋势预测准确率较高.  \n",
    "类似的还有MA, ARMA. 暂未了解.  \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "# 编秩函数\n",
    "def get_rank(data, columns, ascending=True, R=np.zeros(1)):\n",
    "    if data.ndim == 1:\n",
    "        tempdata = np.array(data)[:,None]\n",
    "    if not R.any():\n",
    "        R = np.zeros(tempdata.shape)\n",
    "    for i in columns:\n",
    "        arg = tempdata[:,i].argsort(axis=0)\n",
    "        if not ascending:\n",
    "            arg = arg[::-1]\n",
    "        begin, end = 0, 0\n",
    "        # 找从begin开始相同项, 用end标记最后一个相同项的下一个\n",
    "        while begin < len(arg):\n",
    "            while end < len(arg) and tempdata[arg[end]][i] == tempdata[arg[begin]][i]:\n",
    "                end += 1\n",
    "            for j in range(begin, end):\n",
    "                R[arg[j]][i] = (begin + end + 1) / 2\n",
    "            begin = end\n",
    "    return R.reshape(data.shape)\n",
    "\n",
    "# 作n阶差分折线图\n",
    "def diff_plot(y, n, scale=3):\n",
    "    if y.shape[0] < n:\n",
    "        print('Error')\n",
    "        return\n",
    "    for i in range(n):\n",
    "        y = np.array([y[i]-y[i-1] if i != 0 else y[0] for i in range(y.shape[0])])\n",
    "    plt.figure()\n",
    "    plt.ylim(y.mean()-(y.max()-y.min())*scale, (y.max()-y.min())*scale+y.max())\n",
    "    plt.plot(np.linspace(0, y.shape[0]-1, y.shape[0]), y)\n",
    "    plt.title('{} Order Difference of Data'.format(n))\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "E:\\anaconda3\\lib\\site-packages\\ipykernel_launcher.py:14: RuntimeWarning: divide by zero encountered in double_scalars\n",
      "  \n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "T =  inf\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXIAAAEICAYAAABCnX+uAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3dd3yV9d3/8dcngx1CQhJI2DPsoREVVFBEWXW0zrbWqq1t1VarHdp1a3vfd1tbtVi9tdo6fq21Wq0LcaAVFXGFYQgrREBGAkmYCcgI+fz+OBc2ICPjJCdX8n4+HueRnOtc43My3uc6n/O9rsvcHRERCa+4WBcgIiL1oyAXEQk5BbmISMgpyEVEQk5BLiIScgpyEZGQU5BL1JjZHDP7RiNvc4mZTQi+NzN72My2mtkHwbTvmNkmM6sws86NWVtjMLO2ZvaCmW03s3/Guh6JDQV5M2Rm15lZrpntMbNHajB/dzN7zMw2m9lOM/vAzKY3QqlHq6m3mXkQwBVBGM80s0nV53P3oe4+J7h7CjAJ6O7uY8wsEbgTOMvdO7j75sZ9Fo3iAqAL0NndLzz0QTO71cz2mVl5cCsws3vMLLOmG4jFC7TUjoK8eSoC/ht46FgzmlkqMBfYCwwF0oC7gL+b2QVHWCYheqUec32d3L0DMBKYDTxjZl8/wry9gDXuvjO43wVoAyypY13xdVmukfUCCty98ijzPOHuSUAqcD7QFZhfmzCXJs7ddWumNyJh/sgx5vkVkA/EHTL9x8AngAX3HbgWWAmsDqZNApYD24F7gDeBb1Rbx5XAMmAr8ArQq9pjn1vfIdvvHcyTcMj0HwCbDtQLrAHOBK4CdgP7gQrgcWBnsI4K4N/B/IOIvCBsAVYAF1Vb9yPAfcCsYNkzgdbA74G1wXbvB9oG808A1gM3ASVAMXBFtfW1Be4Ifo7bibxgHlj2JGAesA34CJhwlN/RYGBOMO8S4Jxg+m1EXoD3Bc/xqsMseyvwt0OmxQfb/H1wPwWYCZQGv6uZRN7VAPxP8DPdHWzjnmD6DGAdsAOYD5wa67/3lnyLeQG6NeAvt2ZB/h5w22Gm9wlCMDu470EApgYBlRb8E18AJALfByoJghw4DygMQigB+Bkwr9r6D1rfYbbfm8MHed9g+uDg/hrgzOD7rwNzj7QOoH0QPlcENR0HlAFDg8cfCQJ3HJF3q22APwDPB3UmAS8Avw7mnxA8518GP4OpwC4gJXj83iCAuwXhOZbIC0M3YHMwfxyRF8TNQPphfg6Jwc/xJ0Ar4AygvNrv5VYOCepDlj/s40HN7wffdwa+BLQLnuM/gWerzTuHai/QwbSvBsslEHkh2wi0ifXffEu9qbUiaUT2JA9VXO3xA37t7lvc/VMiIbTU3Z9y931EAm9jtXm/Fcy/zCNv+/8XGGVmvY6wvpoqCr6m1mKZA6YTab087O6V7r4AeJrIi9EBz7n7O+5eBewBvgl8P6izPHgel1Sbfx/wS3ff5+6ziOy1ZptZHJF3JNe7+wZ33+/u89x9D5EQnOXus9y9yt1nA7lEfqaHOgnoAPzG3fe6+7+J7DFfWofnX10Rwc/Q3Te7+9Puvit4jv8DjD/awu7+t2C5Sne/g8gLVHY9a5I6imqvU0KpDDhcrzSz2uMHrKv2fVb1++7uZlb98V7ADDO7o9o0I7I3+slh1ldT3YKvW+qwbC/gRDPbVm1aAvDXaver15ROZC91vpkdmGZE9q4P2OwH96d3EQneNCJ79B8foY4LzewL1aYlAm8cZt4sYF3wwnLAJ/zn51BX3Qh+hmbWjsjnIpOJtFkAksws3t33H25hM7sJ+EZQnwMdOfhFXxqRglxeA75kZrcdEhYXEQm1gmrTqp8qsxjoceCORZKuR7XH1wH/4+6PHWXbdTn15vlE+tEr6rDsOuBNd590lHmq11QGfEqk9bKhltsqI9JX7kekH31oHX9192/WYD1FQA8zi6v2++nJwb+XWgneLXyByO8eIq2RbOBEd99oZqOAhURetOCQ35OZnUrkM5SJwBJ3rzKzrdXml0am1kozZGYJZtaGyJ5jvJm1OcrIkLuI7E39xcy6BvNeCvwU+KG7HylsXwSGmtkXg3V/j8hoiAPuB24xs6FBTclm9rnhcbV4Tl3M7Drgv4BbDnnRqamZwEAzu8zMEoPbCWY2+HAzB9t4ELjLzDKCOrqZ2dnH2lCw7EPAnWaWZWbxZnaymbUG/gZ8wczODqa3MbMJZtb9MKt6n8gHrz8K6p1AJIT/UdsnHyw/mMgHwV2JDM2ESF/8U2BbMIrpvw5ZdBORzyaoNn8lkQ9HE8zsF0T+hiRGFOTN08+I/GPeTKQf+2kw7XM8Mrb6FCJtgKVEPnS7EbjM3Z840gbcvQy4EPhNsMwA4J1qjz8D/Bb4h5ntIDIyZkodnss2M9sJLCbSQ77Q3Y85rPIINZcDZxHpcRcR6en/lkh/90h+TOTDxveC5/EaNe8F/yCo+0MibYzfEhltsw44l8gHmKVE9tB/yGH+H919L3AOkZ9dGfB/wNfcfXkNawC42MwqiIx6eZ7I7+t4dz/wecMfiHyAXUbkw++XD1l+BnBBcKDV3URGIL1E5F3BJ0TeedSlTSZRYkfe4RIRkTDQHrmISMgpyEVEQk5BLiIScgpyEZGQi8k48rS0NO/du3csNi0iElrz588vc/f0Q6fHJMh79+5Nbm5uLDYtIhJaZvbJ4aartSIiEnIKchGRkFOQi4iEnIJcRCTkFOQiIiGnIBcRCTkFuYhIyCnIRURCTkEuIhJyCnIRkZBTkIuIhJyCXEQk5BTkIiIhpyAXEQk5BbmISMgpyEVEQk5BLiIScgpyEZGQU5CLiIScglxEJOQU5CIiIacgFxEJuRoHuZn1MLM3zGyZmS0xs+uD6bea2QYzWxTcpjZcuSIicqiEWsxbCdzk7gvMLAmYb2azg8fucvffR788ERE5lhoHubsXA8XB9+Vmtgzo1lCFiYhIzdSpR25mvYHRwPvBpOvMLM/MHjKzlCMsc7WZ5ZpZbmlpaZ2KFRGRz6t1kJtZB+Bp4AZ33wHcB/QDRhHZY7/jcMu5+wPunuPuOenp6fUoWUREqqtVkJtZIpEQf8zd/wXg7pvcfb+7VwEPAmOiX6aIiBxJbUatGPAXYJm731ltema12c4H8qNXnoiIHEttRq2MAy4DFpvZomDaT4BLzWwU4MAa4FtRrVBERI6qNqNW5gJ2mIdmRa8cERGpLR3ZKSIScgpyEZGQU5CLiIScglxEJOQU5CIiIacgFxEJOQW5iEjIKchFREJOQS4iEnIKchGRkFOQi4iEnIJcRCTkFOQiIiGnIBcRCTkFuYhIyCnIRURCTkEuIhJyCnIRkZBTkIuIhFyNg9zMepjZG2a2zMyWmNn1wfRUM5ttZiuDrykNV66IiByqNnvklcBN7j4YOAm41syGADcDr7v7AOD14L6IiDSSGge5uxe7+4Lg+3JgGdANOBd4NJjtUeC8aBcpIiJHVqceuZn1BkYD7wNd3L0YImEPZESrOBERObZaB7mZdQCeBm5w9x21WO5qM8s1s9zS0tLablZERI6gVkFuZolEQvwxd/9XMHmTmWUGj2cCJYdb1t0fcPccd89JT0+vT80iIlJNbUatGPAXYJm731ntoeeBy4PvLweei155IiJyLAm1mHcccBmw2MwWBdN+AvwGeNLMrgLWAhdGt0QRETmaGge5u88F7AgPT4xOOSIiUls6slNEJOQU5CIiIacgFxEJOQW5iEjIKchFREJOQS4iEnIKchGRkFOQi4iEnIJcRCTkFOQiIiGnIBcRCTkFuYhIyCnIRURCTkEuIhJyCnIRkZBTkIuIhJyCXEQk5BTkIiIhpyAXEQk5BbmISMjVOMjN7CEzKzGz/GrTbjWzDWa2KLhNbZgy/2P3vv0NvQkRkVCpzR75I8Dkw0y/y91HBbdZ0Snr8H41cylTZrzdkJsQEQmdGge5u78FbGnAWo6pV+d2rC7byeqynbEsQ0SkSYlGj/w6M8sLWi8pR5rJzK42s1wzyy0tLa3ThiYMzADgzRUldatURKQZqm+Q3wf0A0YBxcAdR5rR3R9w9xx3z0lPT6/Txnp2bkfftPbMKajbC4GISHNUryB3903uvt/dq4AHgTHRKevIxmen8+7Hm/Whp4hE3a9mLuWKhz+IdRm1Vq8gN7PManfPB/KPNG+0jB+Yzp7KKt5btbmhNyUiLcjL+Rv5y9zVvLGilI9LK2JdTq3UZvjh48C7QLaZrTezq4DbzWyxmeUBpwPfb6A6P3NS3860TojjTbVXRCRKNu3Yzc3/ymNARgcAZuUVx7ii2qnNqJVL3T3T3RPdvbu7/8XdL3P34e4+wt3PcfcGf/ZtEuM5uV9n3lyhIBeR+quqcm568iP27Kvi/suOJ6dXCi8ubqZB3pSMH5jOqrKdrN28K9aliEjIPfTOauYWlvHz6UPol96BqcMzWb6xnMKS8LRXQhnkE7KDYYgFGoYoInW3tGgHt7+8gklDunDpmB4ATB0e+ehvVoj2ykMZ5H3S2tOrczvmqL0iInW0e99+rv/HQpLbJfLbL43AzADomtyGE3qn8GKI+uShDHKItFfmaRiiiNTRr2ctY2VJBXdcOJLU9q0Oemza8ExWbCqnsKQ8RtXVTmiDfEJ2Op/u28+Ha2J61gARCaE3lpfw6LufcOW4Ppw28PMHKE4ZnokZvJi3MQbV1V5og/zkvmm0SojT6BURqZWyij388KmPGNQ1iR9Nzj7sPF06tuGEXqm8uLiokaurm9AGedtW8ZzYJ1WH64tIjbk7P3oqjx27K5lxyWjaJMYfcd5pIzIp2FTByk1Nv70S2iCHSJ+8sKSC9Vs1DFFEju1v733Cv5eXcMuUQWR3TTrqvFOGdY20V0IweiXUQf6fYYjaKxeRoyssKee/X1zG+IHpfH1s72POn9GxDSf0Tg3F6JVQB3m/9PZ0T2mrYYgiclR7KvfzvccX0b51Ar+78D9DDY9l+ohMVpZUUNDE2yuhDnIziwxDLCxjb2VVrMsRkSbqjlcLWFq8g9u/NIKMpDY1Xm7ygfZKE98rD3WQQ6S9snPvfnI1DFFEDuOdwjIeeGsVXzmxJ2cO6VKrZTOS2jCmdyovLi7G3RuowvoLfZCP7deZxHhTn1xEPmfrzr3c9ORH9E1vz8+mDanTOqaPyKSwpIKCTU333CuhD/L2rRM4oXeq+uQichB35yfPLGbzzj3cfclo2rY68lDDozl7WFfimvjoldAHOUSO8lyxqZyibZ/GuhQRaSL+mbuel/I38oOzshnWLbnO68lIasOYPqm8mFfUZNsrzSTII8MQ31J7RUSA1WU7ufWFJZzctzPfPLVvvdc3bUQWH5fuZEUTHb3SLIJ8QEYHMpPbqL0iIuzbX8UNTywiMT6OOy4aSVxczYYaHs3koUF7pYmOXmkWQW5mTMhO553CMvbt1zBEkZbs7tdX8tG6bfzv+cPJ6tQ2KutMT2rNSX07N9nRK80iyAHGD8ygfE8lCz7ZGutSRCRGPli9hXvfKOTC47szbUTmsReohanDM1lVupPlG5tee6U2F19+yMxKzCy/2rRUM5ttZiuDrykNU+axjevfmYQ400m0RGpg3/6qJrlnWR87du/j+08sokdqO/7rnKFRX//kYU23vVKbPfJHgMmHTLsZeN3dBwCvB/djIqlNIsf3SlGfXOQYVmwsZ/ztb3Deve+wpmxnrMuJml88m8/GHbu56+JRdGidEPX1p3Vozcn9OjOrCbZXahzk7v4WcOjhk+cCjwbfPwqcF6W66mRCdgbLinewacfuWJYh0mS9v2ozF9w/j31VzuqynUy7+22eXbgh1mXV27MLN/DsoiKunziA43o2XGNg6vBMVpXtZFlx02qv1LdH3sXdiwGCrxlHmtHMrjazXDPLLS1tmL3mCdmRK33oKE+Rz3tpcTGXPfQBGUmteeaasbx0w2kMyerIDU8s4sYnF1GxpzLWJdbJui27+Pmz+eT0SuGaCf0adFufjV5pYhecaLQPO939AXfPcfec9PTPX1opGgZ1TaJLx9a6apDIIf767hqu+fsChmV15Klvj6V7Sju6dWrL4988iesnDuDZhRv4wh/nsnj99liXWiv7q5wbn1yEA3ddPIqE+IaNtM4dWjO2XxqzFm9sUu2V+j7rTWaWCRB8Lal/SXV34GyIb68spVLDEEVwd37/ygp+/twSJg7K4LFvnERKtQsNJ8TH8f1JA/n7N0/i0737+eJ97/Dnt1dRVdV0Qupo7ptTyIdrtvKr84bSI7Vdo2xz6vBMVpftZGnxjkbZXk3UN8ifBy4Pvr8ceK6e66u3CdkZ7NhdyaJ122JdikhMVe6v4sdP53HPG4VcckIP7v/q8Uc838hJfTvz0vWnMiE7g/9+cRlXPvohZRV7Grni2lm0bht3vbaSc0Zmcd6obo223bOHdiE+zprU6JXaDD98HHgXyDaz9WZ2FfAbYJKZrQQmBfdjalz/NOLjTKNXpEXbtbeSq/86nydz1/O9iQP49ReHH7PtkNK+FQ9cdjy/PHco8z7ezJQZbzN3ZVkjVVw7O/dUcv0/FtK1Yxt+dd6wGl8oIhoi7ZWmdXBQbUatXOrume6e6O7d3f0v7r7Z3Se6+4Dga8xPCp7cNpHjenZiTkFMuzwiMbNl516+/OD7zFlRwn+fN4wbJw2scdCZGV87uTfPXTuO5LaJXPbQ+/zmpeVN7ojp215Ywtotu7jzopEkt01s9O1PG57JJ5t3saSoabRXms2RndVNyM4gf8MOSso1DFFalnVbdnHB/fNYWryD//vK8Xz1pF51Ws/gzI48f924SEvmzY+58P53WbelaVzk/KXFxTyZu55rJvTjxL6dY1LDWUO7RtorTeTUts0yyMcPjIyKebugab4tFGkIS4t28KX75lFWvofHvnEik4d1rdf62rVK4NdfHMG9Xz6Oj0srmDrjbZ7/KDbD7vbtr+KNFSXc9ORH3PjkR4zonswNZw6MSS0Aqe1bRdoreU2jvdIsg3xIZkfSOrTW4frSYsz7uIyL//Qu8XHGU98Zywm9U6O27mkjMpn1vVMZ0KUD33t8IT/850fs2tvwY84r91cxd2UZNz+dxwn/8xpXPPwhry7dyNThmdz/1eNJbOChhscyfUQma7c0jfZK9I9jbQLi4iLDEF9fvon9VU58FE5jKdJUzcwr4sYnPqJX53Y8euWYqJ3xr7oeqe144lsnM+O1ldw7p5D5a7fyx0tHMzSr7hdsOJz9Vc6Ha7YwM6+Il/M3Ulaxl/at4jlzSBemj8jitIFptE6o25V+ou2sIV35yTP5zMwrrteFK6KhWQY5RI7yfHrBej5av61BD9kViaVH3lnNbTOXktMrhQe/lkOndq2OvVAdJcbH8YOzsxnbrzM3PLGI8++dxy1TB/H1sb3rNWqkqspZuG4rL3xUzKzFxZSU76FNYhwTB3Vh+ohMTh+UQZvEphHe1aW0b8W4/mm8uLiIH0/ObtSRM4dqtkF+6oA04gzmrChVkEuz4+7c/soK7pvzMWcN6cLdl45utLAb2z+Nl284jR/+8yNue2Epc1eW8bsLR5LavuYvIu7O4g3bmZlXzMyPiijavptWCXFMGJjO9JFZTByUQfsGOPFVtE0fnsmPns4jf8MOhneP3V550/9J1VGndq0Y1aMTb64o4cZJsftQRCTa9gUH+vxrwQa+fGJPfnXusEZvH6a2b8WfL8/hkXlr+PWs5UyZ8RZ3XTyKsf3SjriMu7OsuJyZeUXMzCtm7ZZdJMYbpw5I5wdnZzNpSBeS2jT+UML6OGtoF37yjDFzcZGCvKFMyM7grtcK2Fyxh84dWse6HJF627mnkmseW8CbBaXcOGkg3z2jf8ze0psZV4zrwwm9U/ne4wv5yp/f59oJ/bnhzAEHHXy0clM5L+QVMzOviFWlO4mPM8b268x1p/fn7KFdSW4XrvCurlO7oL2SV8zNkwfF7HfRrIN8/MB07pxdwNsryzhvdOMdwivSEDZX7OHKRz5k8Ybt/OaLw7lkTM9YlwTAsG7JvPDdU7jthSXc80Yh8z4u4+Ypg3l/1WZm5hWzYlM5cQYn9unMVaf0YfLQrs1qx2raiEx+9FQeeeu3M7JHp5jU0KyDfHi3ZDq3b8WcFSUKcgm1tZt3cfnDH1C07VP+dFkOk4Z0iXVJB2nfOoHbLxjJuP5p/PSZfC7607sAnNA7hdvOGcqU4V3JSGoT4yobxtlDuvLT+MXMWlysIG8IcXHGaQPTebOglKoqj8rVtEUaW/6G7Xz94Q+prKri7988keN7RW+MeLSdO6obx/VM4f3VWxjXvzOZydEfCtnUJLdLZFz/NGbmFXPzlNi0V5rlAUHVjR+Yzpade1m8IVznWRYBmLuyjEseeI/WCXE89e2Tm3SIH9AjtR0XHN+9RYT4AdOGZ7Jh26d8FKPzuTf7ID9tYDoWDEMUCZPnFm3gikc+oHtKW57+zlj6ZyTFuiQ5grOGdCUx3pgVo3OvNPsgT23fihHddTZECY+qKufu11dy/T8WcVzPFJ741sl0TW6e/eXmIrldIqcOSI/ZuVeafZADTBiYzqJ129i6c2+sSxE5qoo9lXz7b/O5c3YB54/uxqNXjonJaVql9qYG7ZVYXNSmRQT5+Ox03OHtQp0NUZquVaUVnHfvO7y+vIRfTB/CnReNbJKHpsvhTRrSJWbtlRYR5CO7dyKlXSJzVqi9Ik3Tv5dv4tx732FzxR7+etUYrjylT0zP3SG1l9w2kdNi1F5pEUEeHxc5DPitYBiiSFPh7tzz75Vc9WguPVPb8cJ3TznqYe7StE0dnknR9t0sbOT2SosIcogMQyyr2NukrnwtLVvFnkq+87cF/P7VAs4dmcVT3x5L95TGuRK8NIwzh3ShVXwcsxr5wswtJshPC64apPaKNAVrynZy/r3v8OrSjfxs2mDuunjUEa9wL+GR3DaR0wamMWtxcaO++49KkJvZGjNbbGaLzCw3GuuMtvSk1gzvlqzx5BJzb6wo4Zx75lJWsYe/XnUi3zi1r/rhzUgs2ivR3CM/3d1HuXtOFNcZVeMHprNg7Va279oX61KkBXJ37n2jkCsf+ZDuKe14/rpTGNdf/fDm5kB75cVGbK+0mNYKRK4aVOUwV8MQpZHt3FPJtX9fwO9eWcEXRmTx9HfG0iNV/fDmqGObRE4bmM5L+Y3XXolWkDvwqpnNN7OrDzeDmV1tZrlmlltaGpv2xqgenejYJkF9cmlUn2zeyRf/bx4v52/kp1MHM+MS9cObu+kjMinevpuF67Y2yvaiFeTj3P04YApwrZmddugM7v6Au+e4e056enqUNls7CfFxnBqcDTEWh9GG1a69lfp51dGbBaV84Y9z2VS+m0evHMM3T1M/vCWYODiDVglxzGyk9kpUTmPr7kXB1xIzewYYA7wVjXVH2/iBkQH7y4rLGZLVMdblNDml5XvIL9rOkg3byd+wgyXF21m35VNyeqVw96WjG+QK7c2Ru3P/m6v43SvLGdgliQe/lqNWSguS1CaR8QPTeWnxRn4+bUiDn0K73kFuZu2BOHcvD74/C/hlvStrIBMODEMsKGnRQe7ubNj2KUuKdkRCu2gH+Ru2U1K+57N5endux4junZg2PIu/vruGaXe/zZ0Xj+L07IzYFR4Cu/ZW8sOn8ngxr5jpIzK5/YIRtGvVrE/9L4cxfUQms5duYsHareT0btjTD0fjr6sL8EzwdjEB+Lu7vxyF9TaIjI5tGJLZkTkrSrlmQv9Yl9MoqqqcNZt3kl+0gyVF21myYQf5RdvZFozeiTMYkJHEKf3TGNotmWFZHRmc1ZGO1S6Ee1FOd655bAFXPPwh35nQj5smDTzouowSsXbzLq7+ay4Fm8q5ZcogrlYrpcWaOLjLZ+2VJh/k7r4KGBmFWhrN+Ox0HnxrFTt27zsorJqDyv1VFJZWkL8hsoe9NAjvnXv3A9AqPo7srklMGdaVIVmR0B7UteMxP3zrm96BZ68dx20vLOW+OR8zf81W7r50tE6vWs1bBaV89/GFADxyxZjPDkKTlqlD6wQmBKNXfjG9YdsrLfL93oSB6dw352PmFZYxeVhmrMups6oq5+PSChau28aiddtYsmE7yzeWs6eyCoC2ifEMyerIBcd3Z2i3ZIZmdWRARhKtEuq2J90mMZ5ff3E4J/ZJ5SfPLGbq3W9z18WjGN/CA8vd+dNbq7j95Ug//IHLcujZWf1wiVyY+dWlm5i/disnNOBeeYsM8uN6pZDUOoE5K0pDFeRbd+5l0bptLFy79bPwLt9dCUBSmwSGZSXztZN7MaxbMkOzkumT1p74BtgLOG90N4Z1S+baxxbw9Yc/4NoJ/bnhzAEtstWya28lP3oqj5l5xUwbnsnvLlQ/XP5j4uAutE6IHBykII+yxPg4xvVPY86KyDDEptjD3Le/ihUbyyOhvXYbC9dtY3XZTiDS0x7UtSNfGJnF6B6dGN0zhb5p7Rv14tL9MyKtllufX8I9bxTy4Zot/PHS0WR0bBmtFndnafEObnryI1ZsKufHkwfx7fHqh8vBOrROYEJ2OrMWN2x7pUUGOUSO8nx5yUYKNlWQ3TX210LcuH03i9YFob12G3kbtrF7X6RFktahNcf17MRFOT0Y3bMTw7sl07517H91bVvF89sLRjCmTyo/ezafqXe/zR8uHs0pA5rnYeebduzmncIy5q4sY25hGSXle+jYJoFHrhjT4ttLcmTTRmTxypJN5H6ylTF9GmavPPZpECPjs/9zNsTGDvLd+/aTv2F7sKcdCe/i7buByIeRQ7t15MtjejG6ZydG9+xEt05tm/Se3peO786I7slc89gCLnvofb57xgCunzigQdo6jaliTyXvr9rM3MIy3ikso2BTBRC5Duy4/mmc0r8zZwzqQnpS6xhXKk3ZxEEZQXulSEEebZnJbRnUNYk5K0r51vh+Dbotd+ftlWW8vmwTC9dtY2nRDiqDczB0T2lLTu/UoEXSiSFZHWmdEL7Dtwd0SeK568bx82eXcPfrK8lds4U/XDKKjKTwtFr27a8ib/023l4ZCe6Fa7dRWeW0TohjTJ9ULji+O+P6pzG4a8dGbWoJnBcAAAq0SURBVGNJuLVvncDp2RnMyt/IL74wtEF2cFpskEPkKM+H3llNxZ5KOjRAq8LdeadwM3fOXsGCtdto1yqekd07cfVpfRndM4VRPTo1q725dq0SuOOikZzUN5WfP5fP1BlzufuSUYxtomf4c4+M+om0Sjbz3qrNVOypxAyGd0vm6tP6ckr/NI7rlaJrZ0q9TBuRyctLNpK7Zgsn9u0c9fW37CDPTudPb61iXmEZZw3tGtV1v7dqM3fOLuCD1VvISm7D/54/nAuO717noX9hcmFOD0Z078Q1j83nq395n+snDuS6M/o3iVZLSflu5hVu/myve+OOSEurV+d2nDMqi1P7p3Fyv850atcqxpVKc3LGoAzaJMbx4uJiBXm05fRKpX2reOYUlEYtyOd/soU7ZxfwTuFmMpJa88tzh3LxCT1C2S6pj+yuSTx/3Sn87Nl87nqtgA+DVktah8Z9B7JzTyUfrN7yWZ97+cZyAFLaJTK2fxqnBDedB0UaUvvWCZwxKINZizfyXw3QXmnRQd4qIY6x/dN4MwrDEBet28Zdswt4s6CUtA6t+Pn0IXzlxJ4t+i15+9YJ3Bm0Wn7x3BKmznibuy8dzUkNsEcCUL57H0uLdhx0KoLC0gr2VzmtEuIY0zuVm6d045T+aQzJVJ9bGtfU4ZnMWryRD9dsifr/QIsOcogMQ5y9dBMfl1bQP6P2o1fyN2znD68V8NqyElLaJXLLlEFcdnIvHRQSMDMuPqEnI7p34trHFvDlB9/jxkkDuWZC/3oF6Zade1lSFDlD44GzNa7ZvOuzxzOSWjOsWzJnDe3CiX06k9NbfW6JrTMGZXDd6f3p1gBnEG3xaTP+s4syl9YqyJdv3MEfZq/k5SUbSW6byA/Pzubysb0b5EPT5mBwZkee/+4p/ORfi/n9qwV8sGYrd100ks7HaLW4O5t27CF/w3aWFP0ntIuC4ZoAPVLbMjQz+aBTEYRptIy0DO1aJfCDs7MbZN0tPnW6p7Sjf0YH5qwo5Run9j3m/IUlFfzhtQJeXFxMh1YJ3HDmAK48pU+zO/lWQ+jQOoEZl4zipL6dufWFJUy9+23+eOlxn42tdXfWbfmU/KLtnwX3kqLtlFXsBcAM+qa154Q+qQzN6siwrGSGZHXUB5PS4rX4IIfISbT+37ufsGtv5RFbImvKdjLj9ZU8t2gDbRPjuXZCf75xah+FSC2ZGV8+sSejenTi2r8v4NIH3+PckVkUbY+cG/3AuWMS4owBXZI4PTsjOHdMRwZndmwSR7SKNDX6rwAmZGfw57mreffjzUwc3OWgx9Zt2cXdr6/kXws3kBhvfPO0vnzrtH6ktleA18eQrI48f904fvZsPq8s2Uj/LkmcMzKLYd2SGZaVzIAuHdTTFqkhBTlwQp8U2ibGM2dF6WdBXrTtU/7470L+mbuOuDjj8pN7850J/ZrVATyxltQmkRmXjG6yJy4TCQsFOdA6IZ6x/Tozp6CETTt2c+8bhfzjg3UAfPnEnlx7en+6tJCz+sWCQlykfhTkgQnZ6by+vIRTf/sGVe5cmNOD685omKFCIiLRpCAPTBrSlfvmfMy4/ml8b+IAHeknIqERlSA3s8nADCAe+LO7/yYa621MXZPbMO+WibEuQ0Sk1up9BicziwfuBaYAQ4BLzWxIfdcrIiI1E41T8Y0BCt19lbvvBf4BnBuF9YqISA1EI8i7Aeuq3V8fTDuImV1tZrlmlltaWhqFzYqICEQnyA83dsw/N8H9AXfPcfec9HRd31BEJFqiEeTrgR7V7ncHiqKwXhERqYFoBPmHwAAz62NmrYBLgOejsF4REamBeg8/dPdKM7sOeIXI8MOH3H1JvSsTEZEaico4cnefBcyKxrpERKR2mv+VgEVEmjkFuYhIyCnIRURCTkEuIhJyCnIRkZBTkIuIhJyCXEQk5BTkIiIhpyAXEQk5BbmISMgpyEVEQk5BLiIScgpyEZGQU5CLiIScglxEJOQU5CIiIacgFxEJOQW5iEjIKchFREKuXkFuZrea2QYzWxTcpkarMBERqZloXHz5Lnf/fRTWIyIidaDWiohIyEUjyK8zszwze8jMUqKwPhERqYVjBrmZvWZm+Ye5nQvcB/QDRgHFwB1HWc/VZpZrZrmlpaVRewIiIi2duXt0VmTWG5jp7sOONW9OTo7n5uZGZbsiIi2Fmc1395xDp9d31EpmtbvnA/n1WZ+IiNRefUet3G5mowAH1gDfqndFIiJSK/UKcne/LFqFiIhI3Wj4oYhIyCnIRURCTkEuIhJyCnIRkZBTkIuIhJyCXEQk5BTkIiIhpyAXEQk5BbmISMgpyEVEQk5BLiIScgpyEZGQU5CLiIScglxEJOQU5CIiIRe1S73VaqNmpcAndVw8DSiLYjmNSbXHRlhrD2vdoNobSi93Tz90YkyCvD7MLPdw16wLA9UeG2GtPax1g2pvbGqtiIiEnIJcRCTkwhjkD8S6gHpQ7bER1trDWjeo9kYVuh65iIgcLIx75CIiUo2CXEQk5EIV5GY22cxWmFmhmd0c63pqysx6mNkbZrbMzJaY2fWxrqk2zCzezBaa2cxY11IbZtbJzJ4ys+XBz/7kWNdUU2b2/eBvJd/MHjezNrGu6UjM7CEzKzGz/GrTUs1stpmtDL6mxLLGIzlC7b8L/mbyzOwZM+sUyxprIjRBbmbxwL3AFGAIcKmZDYltVTVWCdzk7oOBk4BrQ1Q7wPXAslgXUQczgJfdfRAwkpA8BzPrBnwPyHH3YUA8cElsqzqqR4DJh0y7GXjd3QcArwf3m6JH+Hzts4Fh7j4CKABuaeyiais0QQ6MAQrdfZW77wX+AZwb45pqxN2L3X1B8H05kUDpFtuqasbMugPTgD/HupbaMLOOwGnAXwDcfa+7b4ttVbWSALQ1swSgHVAU43qOyN3fArYcMvlc4NHg+0eB8xq1qBo6XO3u/qq7VwZ33wO6N3phtRSmIO8GrKt2fz0hCcPqzKw3MBp4P7aV1NgfgB8BVbEupJb6AqXAw0Fb6M9m1j7WRdWEu28Afg+sBYqB7e7+amyrqrUu7l4MkR0ZICPG9dTVlcBLsS7iWMIU5HaYaaEaO2lmHYCngRvcfUes6zkWM5sOlLj7/FjXUgcJwHHAfe4+GthJ0317f5Cgn3wu0AfIAtqb2VdjW1XLY2Y/JdIWfSzWtRxLmIJ8PdCj2v3uNOG3m4cys0QiIf6Yu/8r1vXU0DjgHDNbQ6SVdYaZ/S22JdXYemC9ux945/MUkWAPgzOB1e5e6u77gH8BY2NcU21tMrNMgOBrSYzrqRUzuxyYDnzFQ3CwTZiC/ENggJn1MbNWRD78eT7GNdWImRmRXu0yd78z1vXUlLvf4u7d3b03kZ/3v909FHuG7r4RWGdm2cGkicDSGJZUG2uBk8ysXfC3M5GQfFBbzfPA5cH3lwPPxbCWWjGzycCPgXPcfVes66mJ0AR58OHDdcArRP6on3T3JbGtqsbGAZcR2aNdFNymxrqoFuC7wGNmlgeMAv43xvXUSPAu4ilgAbCYyP9pkz1s3MweB94Fss1svZldBfwGmGRmK4FJwf0m5wi13wMkAbOD/9X7Y1pkDegQfRGRkAvNHrmIiByeglxEJOQU5CIiIacgFxEJOQW5iEjIKchFREJOQS4iEnL/H6gUBoB3B/ouAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAEICAYAAABcVE8dAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3dd3gU1frA8e9LCaEZuiBd9NJDFUEBkSYoAhYEBQUEsaFerNi4XL0/FVRUlCsXkaYIFqoKgtSogBokFClSpIReQw2Q5Pz+OBNc4iZZssnObvb9PE+e7PR3Z2bfPXvmzBkxxqCUUir3y+N2AEoppQJDE75SSoUJTfhKKRUmNOErpVSY0ISvlFJhQhO+UkqFiVyd8EVkqIh86nYcbhORSiJyUkTyuh2LUp5EpI+I/Bhu23ZLSCd8J4ml/qWIyBmP4Z45sL0mIjJHRI6JyBER+UVE+mb3drKbMWanMaaIMSY5UNsUkQkics45FidEZKWI3HAJyxsRuSonY/SHiFRxYkw937aLyOAAbDd1v55w/taJyOsiEnUJ69guIm1zMs7s4NY+Tmfb+0XkGxFpdwnrCLovlJBO+E4SK2KMKQLsBG71GDc5O7clIs2ARcBS4CqgJPAw0DE7t5PdRCSfi5sf7hybKOBDYHou/JVRzHmPdwIvX0pC8MNwY0xRoDTQF2gK/CQihQOwbTek7uO7gSEi0iHtDDl4nqduux7wPTBDRPrk0LZyXEgnfB9FiMgkpzT0u4g0Tp0gIleIyDQROSgif4rI4xms501gojFmmDHmkLFWGmPu8ljfAyKyxSn9zxaRKzymGRF5REQ2O7G8KiLVRGS5iBwXkS9EJMKZt5WIxIvICyJyyCnZ9PRY1y0isspZbpeIDPWYlloy6SciO4FFHuPyOfP0EZFtThx/pq5bRPKIyEsiskNEDjj7LSrNenuLyE4nrhd9OQDGmBTgM6AEcLlHrPeLyAYROSoi80SksjM+xplltVO66i4iS0XkDmd6cyeWm53htiISl9l6nWk1ROR75xhtEhHP4zdBREaJyLfOvvlZRKr5+B5jgd+B+h7rGywiW511rReR2zym7RCRRs7rXs77qeUM9xeRmT5sM9EY8yvQGVsA6essX01EFonIYec4TRaRYs60T4BKwNfOvn3WGf+liOwTkQQRiRGR2ultV0T6Ovv3hHMePegxLfXcfco5h/aKx69gESnpfDaOi8gvgE/713m/y7H7uI6zLiMij4rIZmCzMy6j4+vPtvcZY94DhgLDRCSPs06vx1hEagKjgWbOfj7mjE/3sxsQxphc8QdsB9qmGTcUSARuBvICrwMrnGl5gJXAECACuBLYBtzkZd2FgGTgxgy23xo4BDQECgDvAzEe0w0wG7gMqA2cBRY6240C1gO9nXlbAUnACGddNwCngOoe0+s67yEa2A90daZVcbY1CSgMFPQYl88Zd9xjXeWA2s7r+4EtTkxFgOnAJ2nW+5GzznrOe6iZzv6YAPzHeZ0XeMjZv3mdcV2dbdV04noJWJZmf13lMfwK8L7z+gVgKzDMY9p7ma3Xee+7sIkxn3OsDnm8/wnAEaCJM30yMDWd93dhnzrDTYHTwG0e83QDrnCOU3fnGJZzpk0CnnJej3Hez8Me0wZltl/TjJ8EfO68vgpohz13SgMxwLuZfFbuB4o6y7wLxGVwrt+CTZaCPTdPAw3TnLuvAPmxn73TQHFn+lTgC+dY1AF2Az9mto+dbV3vrKuNxznyPbYgUdCH45ulbacZf6UzvqYPx7hP2vWTwWc3IHkyUBvK8TeSfsJf4DFcCzjjvL4W2Jlm/ueB8V7WXd45yDUy2P7H2J/aqcNFgPNAFY+T83qP6SuB5zyG3079UHp8aAp7TP8CeDmdbb8LvJPmRL3S28nrnOzHgDuAgmnWsxB4xGO4uvMe8nmso4LH9F+AHunENAH7ZXvM+Z8I9PSYPhfo5zGcB/thruyxvzwTfhtgjfP6O6A/f315LwVuz2y9zgfyhzRx/g/4l0fMYz2m3QxsTOf9pe6PY8AZ5/VbgGRwjsQBXZzX/YDZzusNzvuZ6gzvwEmg6exXbwn/DeD7dJbpCqzK6LOSZv5izvuJ8vGzNxN4wuPcPYNHogQOYL8Q8zrnUw2Paa+RedI9Bhx19tPjHtMN0NpjON3j68e20yb8SNJ8ljM4xn3SW7+3z24g/sKhSmefx+vTQKRTtVEZuELsBdhjzk+uF/CocvBwFEjBlobTcwX2gwqAMeYkcBj7ZZFqv8frM16Gi3hu0xhzymN4h7MNRORaEVkstioqAVt6LpUmnl3egnTW2d1ZZq9TfVHD23twXufj4n2Sdn96xpzWW8aYYtjSV2PgTRFJveZRGXjPY98fwZbiyntfFcuBf4jI5dhqk0lARREphS2Rp1YDZbTeysC1aY55T6BsFt8f2P1eBHgam+zyp04QkftEJM5jW3X46zgtBVqISFlsMvocuF5EqmB/8V2oovJReee9IiJlRGSqiOwWkePAp/z9/LhARPKKyBtO1cRx7BcC6S0jIh1FZIVTbXIM+8XoOe9hY0ySx3DqfiyNPZ88z03P8y09pYwxxY0xNY0xI9NM81xXRsc3q9tOK/X8TN3XGR3jv/Hxs5tjwiHhp2cX8KcxppjHX1FjzM1pZzTGnMYmnDsyWN8e7AkHgNgLaCWxPxuzorhcfBGukrMNsPXhs4GKxpgobF2hpA07vRUbY+YZY9phv8A2Yqtp/vYenG0mcfEX0yUz1jrgJ2x1ANj9/2Ca/V/QGLMsnXWcxv4qegJYZ4w5BywDngS2GmMO+bDeXcDSNNOKGGMe9vP9JRtj3sb+inkEQOx1g4+AgUBJ54tvHc5xMsZswSbCx7FVfyewXzYDsKXCFF+3LyJFgLbAD86o17HHP9oYcxnQi4vPj7Tnxj1AF2cdUdjSLfz9nEJECgDTsL9mLnfe1xxv83pxEHs+VfQYV8mH5TLi+V4yOr7Zte3bsL9YNmV2jPH+GfTls5tjwjnh/wIcF5HnRKSgU8qpIyLXpDP/s0AfEXlGREoCiEg9EZnqTP8M6Csi9Z0PxWvAz8aY7X7E+G8RiRCRFkAn4EtnfFHgiDEmUUSaYD+wPhGRy0Wks/NlchY4ib0+ATAFGCQiVZ0k8hq2XjgpndX5zPkV0Rx70Q3sif586sVBEYkSkW4ei+zH1pd6Wor9cC11hpekGc5svd9gfyXcKyL5nb9rnAts2eEN4FkRicRWnRlsosG5cFknC+8nXSJSQOyF35nYX6HjnUlFscf1mIiUB55Js2jafVsUey4cxl6vei2DzUZg6/kPAknOL7b2vsRrbLPg6cBQESkk9iJ1b1+W9VG6x9ffbTufm4HY6qHnnS/kzI7xfqCCOI0xHFn+7GaHsE34zglwK7Z64E/sxZ2x2BKOt/mXYS/Mtga2icgR7MW2Oc70hcDL2NLPXuxFrR5+hLgP+yHeg714+JAxZqMz7RHgFRE5gb3o/MUlrDcP8JSz3iPYi26PONPGAZ9gq0f+xJZYH/PjPTzrtFA4BczHJqT/ARhjZgDDgKlONcI6Lm7iOhSY6PxUTm1psRT7gYlJZzjD9Tql6PbY47IHu4+HYRNYdvgWe8weMMasx16XWY794NfF/sLxlOn7ScezzrE/gq3aWglc51EF+G/sBcsEJ6bpaZZ/HXjJ2bdPO+vYgf01uh5Ykd6GnX34OPacO4pNWLMzidfTQGz1zj7s9YjxGc59CXw4vlnZ9jHn/F2LrbrqZowZ52wvs2O8CFvA2Sciqb9A/fns+k2cCwcqiIhIK+BTY0wFt2NRSuUeYVvCV0qpcKMJXymlwoRW6SilVJjQEr5SSoUJNzvWuqBUqVKmSpUqboehlFIhZeXKlYeMMaV9nT8oEn6VKlWIjY11OwyllAopInJJdwtrlY5SSoUJTfhKKRUmNOErpVSYCIo6fG/Onz9PfHw8iYmJboeSa0RGRlKhQgXy58+f+cxKqVwnaBN+fHw8RYsWpUqVKogErDO5XMsYw+HDh4mPj6dq1apuh6OUckHQVukkJiZSsmRJTfbZREQoWbKk/mJSKowFbcIHNNlnM92fSoW3oE74SimVm73yCvz2W+C2pwk/A3nz5qV+/frUqVOHW2+9lWPHjmXr+lu1akX16tWJjo6mRo0aDBw40KdtvPZaRs+nUEqFgpkz4V//gulpn1aQgzThZ6BgwYLExcWxbt06SpQowahRo7J9G5MnT2bNmjWsWbOGAgUK0KVLl0yX0YSvVGjbvx8eeAAaNoQhQwK33UwTvoiME5EDIrLOy7SnRcQ4D5JGrJEiskVE1ohIw5wI2g3NmjVj9277eNqTJ0/Spk0bGjZsSN26dZk1axYAw4cPZ+RI+4zlQYMG0bp1awAWLlxIr169Mlx/REQEw4cPZ+fOnaxevRqArl270qhRI2rXrs2YMWMAGDx4MGfOnKF+/fr07Nkz3fmUUsHJGJvsT5yATz6BiIjMl8kuvjTLnAB8gH0M2gUiUhFoB+z0GN0RuNr5uxb40Pnvn3/+E+Li/F7NRerXh3ff9WnW5ORkFi5cSL9+/QDbnn3GjBlcdtllHDp0iKZNm9K5c2datmzJ22+/zeOPP05sbCxnz57l/Pnz/Pjjj7Ro0SLT7eTNm5d69eqxceNG6tWrx7hx4yhRogRnzpzhmmuu4Y477uCNN97ggw8+IM5jf3ibr2TJklnbL0qpHPXxx/D11/DOO1CrVmC3nWkJ3xgTg312ZlrvYB/s7dmhfhdgkrFWAMVEpFy2ROqC1JJ0yZIlOXLkCO3atQNsm/YXXniB6Oho2rZty+7du9m/fz+NGjVi5cqVnDhxggIFCtCsWTNiY2P54YcffEr4qetONXLkSOrVq0fTpk3ZtWsXmzdv9rqMr/Mppdy1dastv7ZuDY8/HvjtZ+nGKxHpDOw2xqxO09SvPLDLYzjeGbfXyzoGAAMAKlWqlPEGfSyJZ7fUOvyEhAQ6derEqFGjePzxx5k8eTIHDx5k5cqV5M+fnypVqpCYmHjh9fjx47nuuuuIjo5m8eLFbN26lZo1a2a6veTkZNauXUvNmjVZsmQJCxYsYPny5RQqVIhWrVp5bUPv63xKKXclJ0Pv3pAvH0yYAHlcuIJ6yZsUkULAi9gnrv9tspdxXh+pZYwZY4xpbIxpXLq0z905uyIqKoqRI0fy1ltvcf78eRISEihTpgz58+dn8eLF7NjxVw+lLVu25K233qJly5a0aNGC0aNHU79+/UzbwJ8/f57nn3+eihUrEh0dTUJCAsWLF6dQoUJs3LiRFStWXJg3f/78nD9/HiDD+ZRSwePNN+Gnn+CDD6BiRXdiyMp3TDWgKrBaRLYDFYDfRKQstkTv+VYqAHv8DTIYNGjQgHr16jF16lR69uxJbGwsjRs3ZvLkydSoUePCfC1atGDv3r00a9aMyy+/nMjIyAyrc3r27El0dDR16tTh1KlTFy4Ad+jQgaSkJKKjo3n55Zdp2rTphWUGDBhAdHQ0PXv2zHA+pVRwiIuzrXG6dQOnrYUrfHqmrYhUAb4xxtTxMm070NgYc0hEbgEGAjdjL9aONMY0yWz9jRs3NmkfgLJhwwafqkHUpdH9qlRgJSZC48Zw5AisXQvZ2Z5CRFYaYxr7Or8vzTKnAMuB6iISLyL9Mph9DrAN2AJ8BDziayBKKZUbvfQS/P67bZ3jduO5TC/aGmPuzmR6FY/XBnjU/7CUUir0LVkCI0bAww9Dx45uR6N32iqlVI5ISLCtcq66yl6wDQZB2x++UkqFsieegPh42zKncGG3o7G0hK+UUtls+nSYOBFefBGCqeGcJnyllMpG+/bBgAHQqBG8/LLb0VxME34GPLtH7tatG6dPn87yupYsWUKnTp28jo+KiqJBgwZUr16dli1b8s033/i0vmXLlmU5HqVU9jMG+vWDU6dsx2jB9vhoTfgZ8OweOSIigtGjR1803RhDSkqK39tp0aIFq1atYtOmTYwcOZKBAweycOHCDJfRhK9U8PnoI5gzB4YNg2C83UUTvo9atGjBli1b2L59OzVr1uSRRx6hYcOG7Nq1i/nz59OsWTMaNmxIt27dOHnyJADfffcdNWrUoHnz5kz38SkH9evXZ8iQIXzwwQcAfP3111x77bU0aNCAtm3bsn//frZv387o0aN55513qF+/Pj/88IPX+ZRSgbNlCzz5JLRpAwMHuh2NdyHRSsfl3pFJSkpi7ty5dOjQAYBNmzYxfvx4/vvf/3Lo0CH+85//sGDBAgoXLsywYcMYMWIEzz77LA888ACLFi3iqquuonv37j7H1rBhQ9502nE1b96cFStWICKMHTuW4cOH8/bbb/PQQw9RpEgRnn76aQCOHj3qdT6lVM5LSoL77rNVOG51jOaLkEj4bkntHhlsCb9fv37s2bOHypUrX+izZsWKFaxfv57rr78egHPnztGsWTM2btxI1apVufrqqwHo1auXzw8n8ezuIj4+nu7du7N3717OnTtH1apVvS7j63xKqew3fDgsXw6TJ0OFCm5Hk76QSPgu9Y58oQ4/rcIejWqNMbRr144pU6ZcNE9cXFymPWSmZ9WqVRf6u3nsscd48skn6dy5M0uWLGHo0KFel/F1PqVU9vrtN/ts2rvugrsz7JfAfUH6wyN0NG3alJ9++oktW7YAcPr0af744w9q1KjBn3/+ydatWwH+9oWQnjVr1vDqq6/y6KO2h4qEhATKly8PwMSJEy/MV7RoUU6cOHFhOL35lFI5JzER7r0XSpeGDz+ELJbxAkYTvp9Kly7NhAkTuPvuu4mOjqZp06Zs3LiRyMhIxowZwy233ELz5s2pXLlyuuv44YcfLjTLfPTRRxk5ciRt2rQBYOjQoXTr1o0WLVpQqlSpC8vceuutzJgx48JF2/TmU0rlnBdegPXrYfx4KFHC7Wgy51P3yDlNu0cOHN2vSmWPRYtsi5xHHoFRo9yJIdu7R1ZKKXWxY8egTx/4xz/sBdtQERIXbZVSKlgYY9vZ79kDy5YFT8dovgjqEn4wVDflJro/lfLfmDG2+eXLL0OTTJ/nF1yCNuFHRkZy+PBhTVLZxBjD4cOHiYyMdDsUpULWL7/A449Dhw72SVahJmirdCpUqEB8fDwHDx50O5RcIzIykgrBfFeIUkHs4EG4804oVw4+/RTy5nU7oksXtAk/f/78ereoUiooJCfbm6oOHLAPNHH72bRZFbQJXymlgsXLL8PChfZB5I0auR1N1gVtHb5SSgWDWbPg9dehf3+4/363o/GPJnyllErH5s22F8xGjeD9992Oxn+ZJnwRGSciB0Rknce4N0Vko4isEZEZIlLMY9rzIrJFRDaJyE05FbhSSuWkU6fg9tshXz6YNg1yQwM3X0r4E4AOacZ9D9QxxkQDfwDPA4hILaAHUNtZ5r8iEoLXspVS4cwY+1za33+HKVMgg66wQkqmCd8YEwMcSTNuvjEmyRlcAaS29esCTDXGnDXG/AlsAULs1gSlVLgbNQo++wxeeQXat3c7muyTHXX49wNzndflgV0e0+KdcX8jIgNEJFZEYrWtvVIqWCxbBoMGQadOtjfM3MSvhC8iLwJJwOTUUV5m83qrrDFmjDGmsTGmcenSpf0JQymlssX+/dCtm63C+eST4H1UYVZluR2+iPQGOgFtzF/9H8QDFT1mqwDsyXp4SikVGElJ0L07HD0Kc+ZAsWKZLxNqsvT9JSIdgOeAzsaY0x6TZgM9RKSAiFQFrgZ+8T9MpZTKWc8/D0uXwujRUK+e29HkjExL+CIyBWgFlBKReOBf2FY5BYDvnee2rjDGPGSM+V1EvgDWY6t6HjXGJOdU8EoplR2++greess+zOS++9yOJucE7ROvlFIqEDZuhGuugTp1bAk/IsLtiHynT7xSSikfnThhb64qWBC+/DK0kn1WaOdpSqmwZAz06webNsH330M49ByuCV8pFZbefdeW6ocNg9at3Y4mMLRKRykVdmJi4Jln4Lbb7P9woQlfKRVW9uyBu+6CatVg/HgQb7eL5lJapaOUChvnz9tkf+KEfaBJVJTbEQWWJnylVNh48kn7iMIpU6B2bbejCTyt0lFKhYV334UPPrBJv0cPt6NxhyZ8pVSuN22aTfS33QbDh7sdjXs04SulcrVly6BXL2jaFCZPhrxh/EgmTfhKqVzrjz+gc2eoWBFmz7Z31IYzTfhKqVzpwAHo2NH2aT93LpQq5XZE7tNWOkqpXOfUKbj1Vti7FxYvtm3ulSZ8pVQuk5wM99wDsbEwfTpce63bEQUPTfhKqVzDGHj8cVtf/8EH0KWL2xEFF63DV0rlGm++Cf/9r+0f59FH3Y4m+GjCV0rlClOnwnPP2efSvvGG29EEJ034SqmQt3Qp9O4NLVvChAm2ZY76O90tSqmQtn49dO0KV14JM2ZAZKTbEQUvTfhKqZC1d69ta1+ggG1rX6KE2xEFN22lo5QKSSdPwi23wOHDtkqnShW3Iwp+mvCVUiEnKcn2a79mjW2C2aiR2xGFhkyrdERknIgcEJF1HuNKiMj3IrLZ+V/cGS8iMlJEtojIGhFpmJPBK6XCjzHw8MO2CufDD+Hmm92OKHT4Uoc/AeiQZtxgYKEx5mpgoTMM0BG42vkbAHyYPWEqpZT1f/8HY8fCiy/CAw+4HU1oyTThG2NigCNpRncBJjqvJwJdPcZPMtYKoJiIlMuuYJVS4W3SJHj5Zbj3Xnj1VbejCT1ZbaVzuTFmL4Dzv4wzvjywy2O+eGfc34jIABGJFZHYgwcPZjEMpVS4WLgQ+vWD1q1tCT+cHj6eXbK7Waa3Q2C8zWiMGWOMaWyMaVy6dOlsDkMplZusXQu33w41atgO0SIi3I4oNGU14e9Prapx/h9wxscDFT3mqwDsyXp4Sqlwt3EjtGsHRYrAnDkQFeV2RKErqwl/NtDbed0bmOUx/j6ntU5TICG16kcppS7VH3/YKhyABQvsk6tU1mXaDl9EpgCtgFIiEg/8C3gD+EJE+gE7gW7O7HOAm4EtwGmgbw7ErJQKA5s3w4032jb3S5ZAzZpuRxT6Mk34xpi705nUxsu8BtBOSZVSftm61Sb7c+fsE6tq1XI7otxB77RVSgWVbdtssk9MhEWLoE4dtyPKPTThK6WCxvbtNtmfOmWbYUZHux1R7qIJXykVFHbssMn++HGb7OvXdzui3EcTvlLKdbt22WR/9KhN9g21F64coQlfKeWq+Hib7A8ftk0vtefLnKMPQFFKuWbPHtvO/uBBmD8frrnG7YhyNy3hK6VcsXevLdnv3WuT/bXXuh1R7qcJXykVcPv22ZL97t0wbx40a+Z2ROFBq3SUUgF14AC0aWMv1M6dC9df73ZE4UNL+EqpgDl40Jbst2+3HaG1aOF2ROFFE75SKiAOHbIl+23b4Ntv4YYb3I4o/GjCV0rluMOHoW1b2yHaN9/Yi7Uq8DThK6Vy1JEjtj/7jRth9mxbylfu0ISvlMoxR49C+/bw++8wa5Z9rdyjrXSUUjkitZ392rUwYwZ06OB2REpL+EqpbLdpE9x0k627//prLdkHC034Sqls9fPPcMstkDevfVKV9o0TPLRKRymVbb791lbjFCsGy5Zpsg82mvCVUtli3Djo0sU+jvCnn6BaNcAYW6+jgoJW6Sil/GIMvPYavPQStL/uBF/1+Zaib8bCqlUQFwcFC9o+kJXrNOErpS7dqVOwZg3JK+N44sMajFp/Iz3zTGHcst5ELDsPkZFQty7ccQc0aAApKZBHKxTcpglfKZWxAwdsSX3Vqr9K7X/8QaKJoBefMo0beabS57xx52/kaTDOJvjq1SGfppdg49cREZFBQH/AAGuBvkA5YCpQAvgNuNcYc87POJVSgRQbC0OH2gS/Z89f4ytXhvr1Oda1D12+6U/M76UY8bZh0JPdge5uRat8lOWELyLlgceBWsaYMyLyBdADuBl4xxgzVURGA/2AD7MlWqVUztu1y7arFLF9IjRoYJ8oXr8+lCjB7t32JqpNf8Bnn8Hdd4vbESsf+fubKx9QUETOA4WAvUBr4B5n+kRgKJrwlQoNiYm23v30adugvlatiyZv2GBvqDp61HZv3LatS3GqLMnyVRRjzG7gLWAnNtEnACuBY8aYJGe2eKC8t+VFZICIxIpI7MGDB7MahlIquxgDDz8Mv/4Kn3zyt2S/bJl9WMm5cxATo8k+FGU54YtIcaALUBW4AigMdPQyq/G2vDFmjDGmsTGmcenSpbMahlIqu4waBRMmwMsvQ9euF036+mub4EuWtIm/QQN3QlT+8aedVFvgT2PMQWPMeWA6cB1QTERSq4oqAHvSW4FSKkjExMCgQdCpk71Y62HsWJv/69SxN1RdeaU7ISr/+ZPwdwJNRaSQiAjQBlgPLAbudObpDczyL0SlVI7atQu6dbOZ/NNPL7SXNwZefRUeeMB2frZoEZQp43Ksyi/+1OH/DHyFbXq51lnXGOA54EkR2QKUBD7OhjiVUjkh9SLtmTMwcyZERQGQnAyPPAJDhsB999kHlxQp4nKsym9+tdIxxvwL+Fea0duAJv6sVykVAJ4XaWfMgJo1AdsCp2dPmDsXBg+23SaItrzMFfRWOKXClZeLtOvW2Zc7d8Lo0fDgg+6GqLKXJnylwpGXi7RffQV9+kDRorB4sW2CqXIX7c1IqXCzaxfceeeFi7TJJg/PP2+v29atCytXarLPrbSEr1Q4OXMGbr/dXqydOZMjyVHccwvMm2db47z/PhQo4HaQKqdowlcqXKRepI2NhZkzWZtUk67X2AL///4HAwa4HaDKaZrwlQoXH3wAEyfCkCF8cbYLfZvaVphLl0KzZm4HpwJB6/CVCgdLl8KgQSR36sJzZ4bSvbvt/HLlSk324UQTvlK5nXMn7ZGqjbj59JcMf1N46CHbEqdcObeDU4GkVTpK5WbORdrVp67itgJL2P1jfj76CPr3dzsw5QZN+ErlVs5F2qmx1bg/4lOKp+Rj6VJo2tTtwJRbtEpHqVwq6b1RPDOxNnczlYbX5GPlSk324U5L+ErlQodn/0SPQTVYQFsefsjw7ntCRITbUSm3acJXKpeJm7uX226rwB4px9iRZ+g3sKDbIakgoVU6So3nKW0AABZDSURBVOUSxsDHH57juluKc97k44epezTZq4towlcqFzh8GO7seJL+j0TQ1Cwjdvw6mtxVxe2wVJDRhK9UiFu4EKKvPsPX8yIYHjmEBdNPULb3TW6HpYKQJnylQtTZs/D0oCTatoWiR3ewonZ/ntlwP3lu6+J2aCpI6UVbpULQ+vVwz51nWb2hAA/xIW8/toNCb41Fm+KojGgJX6kQYox9UFWj+sns3niC2YXv5sOZV1Bo5Bua7FWmtISvVIjYvx/u75PMnO/y0oH5jG/wPmWn/xeqVHE7NBUiNOErFQK+/Rb63pfM8aNJjOSfDBwUgbwxU0v16pJolY5SQezMGRg40D56ttyx9cQWbc1js9ohI97WZK8umV8lfBEpBowF6gAGuB/YBHwOVAG2A3cZY476FaVSYSguDu65O4UNG/MwiBG81mgmkV9+BpUrux2aClH+lvDfA74zxtQA6gEbgMHAQmPM1cBCZ1gp5aOUFHjrLWjSxHBsy2Hm044RT+0h8scFmuyVX7Kc8EXkMqAl8DGAMeacMeYY0AWY6Mw2Eejqb5BKhYvdu6F9e3jmGehkvmZN4Wa0m/WY/QbQKhzlJ3+qdK4EDgLjRaQesBJ4ArjcGLMXwBizV0TKeFtYRAYAAwAqVarkRxhK5Q7TpsEDDxjOnjjHRzxKv4brkC8WaqleZRt/qnTyAQ2BD40xDYBTXEL1jTFmjDGmsTGmcenSpf0IQ6nQtm8f9OgBd94J1c6uZ1VSXfo/VQz5IUaTvcpW/iT8eCDeGPOzM/wV9gtgv4iUA3D+H/AvRKVyp5QUGDMGatY0zJiWzL8j/o9lETfyj1lvaRWOyhFZTvjGmH3ALhGp7oxqA6wHZgO9nXG9gVl+RahULrRhA9xwAzz4INRLiWNNUi2GNPue/HG/QufOboencil/b7x6DJgsIhHANqAv9kvkCxHpB+wEuvm5DaVyjcREeP11eP11Q5G8Z/g47xP0zTMNGfc29OkDIm6HqHIxvxK+MSYOaOxlUht/1qtUbrRkiS3R//EH9Cw6mxEnHqBMz/YwYiOU8dq2QalspV0rKJXDDh+2zSzHj4eqRQ4wj160L70Vpk2Gdu3cDk+FEe1aQakcYgxMnmwvyk6amMJzBUey7sxVtB/cCNau1WSvAk5L+ErlgG3b4OGHYf58aHLZRr5P6UG9eoVgzE9Qt67b4akwpSV8pbLR+fMwbBjUqWNYtuQs7+cbxDJzHfVGPQg/abJX7tISvlLZ5OefYcAAWLMGuhZdxPtnelPhjqYw8ne44gq3w1NKE75S/jp+HF58EUaNMlxR6Bgz6EvXYr/B5A/h1lvdDk+pC7RKR6ksSkmBTz+FWrUMo0YZHi04jvWnq9L1n1XtQ2c12asgoyV8pbIgJgaeegpiY6Fh1FammZ5cW/08fLQQGjVyOzylvNISvlKXYPNmuP122y3Cvt1JTCr7LL+eqcu1b3eHX37RZK+CmpbwlfLBkSPwyiswahQUKACvPryHJ6ddT6HzCbBgPrRo4XaISmVKE75SGTh3zib5V1+FhATo1w9euXExZQd0hhIlYMlPULOm22Eq5ROt0lHKC2Ng+nSoXRuefBKuucY+Y3bMdRMoe197qFYNli/XZK9CiiZ8pdL49VdbR3/HHbb6Zu5cmPedoe7s/4O+fe3EmBhtW69CjiZ8pRw7d0KvXtCkCWzaBKNH21J9h7ZJtp+El16yM8yZA5dd5na4Sl0yrcNXYe/ECXjjDRgxwg4//zwMHuzk9NOn7fMHv/7ajnztNe2zXoUsTfgqbCUlwccfw5AhcOAA9Oxp83mlSs4MBw/am6d++QU++AAefdTVeJXylyZ8FXZSUmDaNPj3v+H336F5c/jmG3th9oKtW6FjR9i1y858222uxatUdtE6fBU2UlLgiy8gOhruuguSk20uj4lJk+xjY+G66+yTSxYu1GSvcg1N+CrXS06Gzz+3PRN3724T/5QpsG6dvWv2oir5uXOhVSsoVAiWLbOJX6lcQhO+yrWSk21ir1vXXncFmDrVPmyqRw/ImzfNAuPG2Tr7f/zDtrGvXj3gMSuVkzThq1wnORk++wzq1IF77oE8eWwJf+1aW8L/W6I3xlbo9+sHbdrA0qVQtqwrsSuVkzThq1wjKcl2V1y7tm1xky8ffPmlfSDJXXfZxO91oQEDYOhQ6N3bXr0tWjTQoSsVEH4nfBHJKyKrROQbZ7iqiPwsIptF5HMRifA/TKXSl5QEn3xiE/2990JEBHz1FaxeDXfemU6iBzh1Crp2hbFj7U1V48dD/vwBjV2pQMqOEv4TwAaP4WHAO8aYq4GjQL9s2IZSf5OUBJMmQa1acN99ULCgbXUTF2e7RUg30YMt9t9wg71IO3q07R1Nb6hSuZxfCV9EKgC3AGOdYQFaA185s0wEuvqzDaXSSkqCCRNsv2W9e0PhwjBjBvz2m211k2Gi37/fVuE0aAB//gkzZ8KDDwYqdKVc5W8J/13gWSDFGS4JHDPGJDnD8UB5bwuKyAARiRWR2IMHD/oZhgoHJ07AyJG2EU3fvraqfeZMm+i7ds0k0ScmwrBhcPXVturmiSdgyxZ9DKEKK1lO+CLSCThgjFnpOdrLrMbb8saYMcaYxsaYxqVLl85qGCoM7NgBTz8NFSrYPF2uHMyaBStXQpcumdTEGGMr9GvVsn3htGplb68dMQKKFw/UW1AqKPjTtcL1QGcRuRmIBC7DlviLiUg+p5RfAdjjf5gqHC1fDu+8Y/ulB+jWDQYNsr1Z+mTlSrvADz/Yxvjffw9t2+ZYvEoFuyyX8I0xzxtjKhhjqgA9gEXGmJ7AYuBOZ7bewCy/o1RhIynJtplv2tTe5Dp/vn0AybZt9iYqn5L9nj22zueaa2DjRvjf/2DVKk32KuzlROdpzwFTReQ/wCrg4xzYhspljh2zrSPff9/2S3/VVfZ1nz5QpIiPKzlzBt5+2/Z1fP48PPMMvPACREXlZOhKhYxsSfjGmCXAEuf1NsDXH90qzG3dCu+9Z3s1OHXKtpQcORI6dfJyR2x6jLF9JgwebL8t7rjDXqCtVi1HY1cq1Gj3yCrgjLE9VL7zDsyebe+I7dHDVrc3aHCJK/v5Z7vg8uV24UmT7LeGUupvNOGrgDl3znZPPGKErVIvWdLWuDzySBYeD7trl3001eTJtt+bcePs3Vc+/yxQKvxowlc5bssW+2SpCRNg3z6oUcNeR+3Vy/ZC7LOUFNvyZto0W++TkgIvvgjPPaf93yjlA034KkckJtrmlGPHwuLF9qaoW26xzwK/6aZMbpLydOiQbaozdy7Mm2cfOyhie0MbNgwqV87R96FUbqIJX2WrtWttkv/kEzh6FKpWhf/8x7a2Ke/1nus0kpPtE6fmzrV/v/5qK/1LlbLfFB07Qvv2oDfrKXXJNOErv508aRvJjB1rr6FGRNinAvbvD61b+1CaP3Dg4lL84cO2FH/ttbbb4g4doFEjrZ9Xyk+a8FWWGGML3x99ZJP9yZO2M7MRI2wXxaVKZbBwcjL88stfpfiVK+0Ky5Sx9T4dOthSfMmSAXs/SoUDTfjqkhw5YhvGfPSRrb4pVMg+Rap/f2jWLJ1+bc6etf3X/PYbLFhgS/NHj9qif9Om8MortqqmQYNLqNxXSl0qTfgqU8bYp/6NHWv7ITt71tawjB4Nd98Nl13mMfOxY7ZD+rg42/YyLg7Wr7d9JoBtQtmli03wbdtCiRKuvCelwpEmfJWu33+3/ddMmWL7somKso997d8fGtQ3sHs3LF11cXL/88+/VlCunC21d+oE9evbv2rVtBSvlEs04auLbN9u6+SnTLEPhcqTB9q0TuFf/fdyZ9kfKbQ+Fp51Evzhw3YhEdvPfJMmfz1cpH59uPxyV9+LUupimvAV+7ef4cuPjzNlegTL1ts+4puV3szIanO46+wnXL44DhYk25kLFLBdDd9221+JPTr6Eno4U0q5RRN+bmeMrVffseOiv4Sth5ix5kqm7G7JgnMtSeFy6rKG1xhOjzxfUrXgebiiMlSuCZU7QPXqNsFXr64P+lYqRGnCz40SE+1DXseNsw3jT5wA4AyRfMstTMnTi29NR86aAlQtepDBzVZw901HqNO8GFR+GK541fZoppTKVfRTnZusXm07rfn0U9vssUoVzvfqy8KkG5jyRyNmxFbgxKm8XF7a8GB34Z57oEmT0ojoXatKhQNN+KEuIcFeYR071t7AFBFBYpfuLKo3iFk76jP9S+HQIdvCplt324zyxhtFb1pVKgxpwg9FqR3Kf/yxbRh/5gxHa13Pt72+Y9bxVnw3twAnv7TXUW++Ge65x968WqCA24ErpdykCT+U7N0LEyfaRL9lC9uL1GFWgwnMOteBmFVFSV4vlC1rE3yXLrYfm8hIt4NWSgULTfjBLikJ5syBsWMx385hVUo0syq9yKxKnVm9swQsg1q14NlnbZK/5hq9r0kp5Z0m/GC1eTN8/DHnJnzG0v3VmVXobmYXmcSu48XIEw/XXQdvPmaT/NVXux2sUioUaMIPJps3w4wZJHwxj+9WlmKWdGVO3pdIoAgFjaH9jcK/u9ieCrQ7eKXUpcpywheRisAkoCyQAowxxrwnIiWAz4EqwHbgLmPMUf9DzYWMgdWrSZ42k9jP/mDetquYx038zJMkk49SxZO5vUteunSBdu3k0h4HqJRSafhTwk8CnjLG/CYiRYGVIvI90AdYaIx5Q0QGA4OB5/wPNZdISYHly4mfuJD5M04x71BDFvAYRyiJiKFR3XMMvjUfHTpAs2Z5tfmkUirbZDnhG2P2Anud1ydEZANQHugCtHJmmwgsIdwT/rlznJkXQ8zo9cxbUoB5p5uzniEAlIs6ReeOeWnf2ZbiS5XStpNKqZyRLXX4IlIFaAD8DFzufBlgjNkrImWyYxuhxpw8xfqxy5j3yQHmrSlLTNL1JNKWAnnO0bLuEfrelchNXSKpU6ew94eGKKVUNvM74YtIEWAa8E9jzHHxMXuJyABgAEClSpX8DSMoHFx/kEWjNjDvm/PM31mD3bQDoGbUbh66YS833V+elu0KUKhQWZcjVUqFI78Svojkxyb7ycaY6c7o/SJSzindlwMOeFvWGDMGGAPQuHFj408cbtm1YjcxE/8kZkkKMdvKs/FcNaA0xeUYbatt46aux2n/6NVUrFre7VCVUsqvVjoCfAxsMMaM8Jg0G+gNvOH8n+VXhEHCpBg2z/+TmM/iifkxDzE7q7AjuQJQnigSaF5mE30b7eKGO0vT+N6a5M3f0O2QlVLqIv6U8K8H7gXWikicM+4FbKL/QkT6ATuBbv6F6I7kc8msm7GZmC/3EbOiADF7ruKAuRK4kjJykJZXbOGppltpeVdZ6nS9irwRTdwOWSmlMuRPK50fgfQq7Ntkdb1uOXfyHL9N/YOY6YeIWVmIHw9UJ4EaQA0q543npqqbaNl8Iy3vqcDV7aogefTOJ6VUaMlVd9qePX6WhF3HSdhzimN7TpOwP5GEg2c5djCJhKPJHDsKCceFhJN5OXYqPwmJESScjeTY+cLsSyrFGeoAUCNiK91rrKbljXlpcW8VKjWrAFRw980ppZSfQjrhz33lVwb9X0mOJRUhIaUoiRQESjt/fyekcBkniMp7gmL5TxEVkUj5oieoXfgwZUps4/rWkTTvXY0ytasB1QL5VpRSKseFdMIvfkVBosvsp1iReKKKpFCsGEQVE4qVykdUqfxElSlAsXIFiSpXiGIVi1L0iqLkyRcFRLkdulJKBVxIJ/ym/evwRX+3o1BKqdCgPacrpVSY0ISvlFJhQhO+UkqFCU34SikVJjThK6VUmNCEr5RSYUITvlJKhQlN+EopFSbEGPe7oheRg8COLC5eCjiUjeEEgsYcGKEWc6jFCxpzoKQXc2VjjM89OQZFwveHiMQaYxq7Hcel0JgDI9RiDrV4QWMOlOyKWat0lFIqTGjCV0qpMJEbEv4YtwPIAo05MEIt5lCLFzTmQMmWmEO+Dl8ppZRvckMJXymllA804SulVJgImYQvIh1EZJOIbBGRwV6mFxCRz53pP4tIlcBHeVE8FUVksYhsEJHfReQJL/O0EpEEEYlz/oa4EWuamLaLyFonnlgv00VERjr7eY2INHQjTo94qnvsvzgROS4i/0wzj+v7WUTGicgBEVnnMa6EiHwvIpud/8XTWba3M89mEentYrxvishG57jPEJFi6Syb4TkU4JiHishuj2N/czrLZphfAhzz5x7xbheRuHSWvfT9bIwJ+j8gL7AVuBKIAFYDtdLM8wgw2nndA/jc5ZjLAQ2d10WBP7zE3Ar4xu39myam7UCpDKbfDMwFBGgK/Ox2zGnOk33Ym1GCaj8DLYGGwDqPccOBwc7rwcAwL8uVALY5/4s7r4u7FG97IJ/zepi3eH05hwIc81DgaR/OmwzzSyBjTjP9bWBIdu3nUCnhNwG2GGO2GWPOAVOBLmnm6QJMdF5/BbQREQlgjBcxxuw1xvzmvD4BbADKuxVPNuoCTDLWCqCYiJRzOyhHG2CrMSard23nGGNMDHAkzWjPc3Yi0NXLojcB3xtjjhhjjgLfAx1yLFCHt3iNMfONMUnO4AqgQk7HcSnS2ce+8CW/5IiMYnby113AlOzaXqgk/PLALo/heP6ePC/M45yUCUDJgESXCad6qQHws5fJzURktYjMFZHaAQ3MOwPMF5GVIjLAy3RfjoVbepD+hyPY9jPA5caYvWALCEAZL/ME6/6+H/tLz5vMzqFAG+hUQ41Lp9osWPdxC2C/MWZzOtMveT+HSsL3VlJP257Ul3kCTkSKANOAfxpjjqeZ/Bu2+qEe8D4wM9DxeXG9MaYh0BF4VERappkerPs5AugMfOllcjDuZ18F3f4WkReBJGByOrNkdg4F0odANaA+sBdbRZJW0O1jx91kXLq/5P0cKgk/HqjoMVwB2JPePCKSD4giaz/vso2I5Mcm+8nGmOlppxtjjhtjTjqv5wD5RaRUgMNMG9Me5/8BYAb2564nX46FGzoCvxlj9qedEIz72bE/tTrM+X/AyzxBtb+di8adgJ7GqUhOy4dzKGCMMfuNMcnGmBTgo3RiCap9DBdy2O3A5+nNk5X9HCoJ/1fgahGp6pTkegCz08wzG0htwXAnsCi9EzIQnPq3j4ENxpgR6cxTNvU6g4g0wR6Pw4GL8m/xFBaRoqmvsRfp1qWZbTZwn9NapymQkFot4bJ0S0PBtp89eJ6zvYFZXuaZB7QXkeJOdUR7Z1zAiUgH4DmgszHmdDrz+HIOBUya60u3pROLL/kl0NoCG40x8d4mZnk/B+JKdDZdzb4Z29JlK/CiM+4V7MkHEIn9Ob8F+AW40uV4m2N/Fq4B4py/m4GHgIeceQYCv2NbBawArnM55iudWFY7caXuZ8+YBRjlHIe1QOMgODcKYRN4lMe4oNrP2C+jvcB5bImyH/Ya00Jgs/O/hDNvY2Csx7L3O+f1FqCvi/FuwdZ1p57Pqa3irgDmZHQOuRjzJ855ugabxMuljdkZ/lt+cStmZ/yE1PPXY16/97N2raCUUmEiVKp0lFJK+UkTvlJKhQlN+EopFSY04SulVJjQhK+UUmFCE75SSoUJTfhKKRUm/h+lUsisCcxQKQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# 15.5.2 税收收入AR预测模型\n",
    "\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "at = np.array([15.2, 15.9, 18.7, 22.4, 26.9, 28.3, 30.5, 33.8, 40.4, 50.7, 58, 66.7, 81.2, 83.4])\n",
    "\n",
    "# 1.检验序列是否平稳\n",
    "# 计算秩\n",
    "Rt = get_rank(at, columns=[0], ascending=True)\n",
    "# 计算Spearman相关系数qs\n",
    "qs = 1 - 6*(sum([(i+1-Rt[i])**2 for i in range(Rt.shape[0])]))/(at.shape[0]*(at.shape[0]**2-1))\n",
    "# 计算统计量T\n",
    "T = qs*np.sqrt(at.shape[0]-2) / (np.sqrt(1-qs**2))\n",
    "print('T = ', T)\n",
    "# 按照书上的一些判断可以知道序列非平稳.\n",
    "\n",
    "# 先画折线图, 书上说可由此看出一阶差分时间序列是平稳的. 不太理解.\n",
    "diff_plot(at, 1, scale=1)\n",
    "\n",
    "# 用差分方程模型(似乎就是所谓AR模型?)\n",
    "y = at[2:]\n",
    "x = np.array([at[1:-1], at[0:-2], np.ones(at.shape[0]-2)]).T\n",
    "C = np.matmul(np.linalg.inv(np.matmul(x.T, x)), np.matmul(x.T, y))\n",
    "\n",
    "pred_at = at[:3]\n",
    "for i in range(15):\n",
    "    pred_at = np.r_[pred_at, [C[0] * pred_at[pred_at.shape[0]-1] + C[1] * pred_at[pred_at.shape[0]-2] + C[2]]]\n",
    "\n",
    "plt.figure()\n",
    "plt.plot(np.linspace(0, at.shape[0]-1, at.shape[0]), at, c='red', label='Raw Data')\n",
    "plt.plot(np.linspace(0, pred_at.shape[0]-1, pred_at.shape[0]), pred_at, c='blue', label='Pred Data')\n",
    "plt.legend()\n",
    "plt.title('The Comparison Between Raw Data and Pred Data')\n",
    "plt.show()\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
